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Abstract 


The aim of the present work is to develop a computer program to compute the 
torque and power generated by a Pelton Wheel with a given bucket surface con- 
figuration. The water jet striking the Pelton Wheel buckets is modelled as a series 
of particles. The force and the torque exerted by every particle, both on impact 
and subsequently during the motion on vane, are computed, using Newtonian me- 
chanics, at every instant of time and integrated over entire time span until they 
lose contact from the bucket. The time history of all particles are then used to 
obtain the total torque obtained by water jet, as a function of time. The bucket 
shape is taken to be arbitrarily specified in local spherical coordinates r = / (6, (f>) 
for which data are supplied as rectangular grid (M x N) in 6 and (f> having the 
value of r at every grid point. The aim is to help in the optimum design of Pelton 
Wheel buckets. 



Chapter 1 
Introduction 


The purpose of a hydraulic turbine is to convert kinetic and potential energy of 
water into mechanical work. On the basis of the difference between the lower and 
higher water levels, known as head, three principal types of turbines - impulse 
turbine, the radial or mixed-flow reaction turbine, and the axial flow, or propeller, 
reaction turbine can be used. 



Figure 1.1: Schemetic diagram of Pelton wheel turbine 
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Among the impulse water turbines, the Pelton turbine as shown in figure 1.1 is 
the only type being mostly used these days. It is also called the free jet turbine and 
operates under a high head and, therefore, requires a comparatively less flow rate 
of water. The head of the supply water is converted into one or more high-velocity 
free jets which are directed against buckets mounted on the rim of a wheel. The 
impact of water on the surface of the bucket produces a force which causes the 
wheel to rotate, thus supplying the mechanical power to the shaft. 


r Y 
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Figure 1.2: schemetic diagram of pelton wheel bucket 

In general, each bucket consists of two hemispherical cups separated by a 
sharp edge at the centre (see figure 1.2 ). The water jet impinges at the centre 
of the bucket and is divided by the sharp edge, without shock, into two parts 
moving sideways in opposite directions. The water after impinging on the buckets 
is deflected through an angle of about 165° instead of 180°as shown, so that it may 
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not strike the back of the other bucket and retard the motion of the wheel. The 
change in the momentum due to this change of directions transmits an ’impulse’ 
torque on the bucket which translates to work output. 

The efficiency of the wheel is given as 

work done on the bucket 

TJ = 

energy supplied to the bucket 

The major problem is to find the power developed by a particular machine. 
The shape of the bucket plays a very important role in creating the torque. If 
the bucket surface is not proper, the water jet is not deflected smoothly with the 
result that there is needless dissipation of energy due to turbulence, and other 
irregularities in the bucket. More important, if the bucket is improperly shaped, 
the water flow does not undergo a sufficient deflection and hence generates less 
torque. So the optimum surface which gives maximum efficiency is required. 

The efficiency of the Pelton wheel is generally very high (>90%). Therefore 
if improvement has to be made, we need to have analytical methods that can de- 
tect the effects of subtle changes in the bucket shape. In most text books, two 
dimensional analysis is presented using the velocity traingle for the calculation of 
work done on the bucket. This simple analysis using the velocity triangle is not 
adequate since it does not have any consideration of surface geometry except the 
deflected angle. Through experiments we can also analyse the effect of the bucket 
surface on the power developed by the wheel. But changing bucket parameters 
and doing experiment requires a lot of time as well as resources. 

It would therefore, be of great help to the designer if a computational method 
can be found to theoritically evaluate the efficiency of a proposed bucket shape. 
The speed and relatively low cost of computational solutions will allow the de- 
signer to evaluate many possible bucket shapes before actually constructing pro- 
totypes for experimental testing. This problem, though an important one, has no 
published work. The commercial importance of the problem has ensured that any 
research on this problem is kept secret. 
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The general problem can be stated deceptively simply: Given a particular 
bucket shape, can its performance in terms of efficiency, torque, force and power 
generation can be predicted? The problem is a very difficult one. Any attempt to 
obtain answers through the traditional routes of fluid mechanics seems to be cur- 
rently impossible- for the jet/bucket interaction involves a highly unsteady, even 
turbulent, flow in a complex geometry involving free surfaces. To attempt to solve 
for flow using Euler or Navier Stokes equations seems to be completely outside 
the scope of present-day computational fluid dynamics. 

In this thesis a method, far simpler and probably more likely to suceed, is at- 
tempted. The jet is treated as a steady stream of perfectly plastic (i.e, deformable 
) “particles” which strike the bucket and slide on its surface imparting momentum 
and energy to the bucket. The motion of each such “particle” is taken to be inde- 
pendent of others and is computed using the Newtonian laws governing particle 
motion. These are coupled ordinary differential equations ( much more numeri- 
cally tractable than the partial differential equations of fluid flow ) which can be 
solved by the Runge- Kutta method. However, considerable difficulty is intro- 
duced into the problem due to the generally complex bucket geometry involved. 
Once the time history of each particle is known from the solution of the dynamical 
equations, the effects of all particles can be integrated to obtain the total picture 
regarding the efficiency, torque etc., generated by the bucket. 

In the present work, we have developed a computer progarm to analyse the 
performance of the buckets of general shapes of 2 and 3 dimensions. The torque 
imparted to the bucket is computed at every instant of time and integrated over 
entire time span. 

The geometry of the bucket is usually very' complex and is very important 
in the calculation of torque. In the present work, the surface of the bucket is 
represented in spherical coordinates (r, 6, <?) having uniform grid in 6 and <f>, with 
the value of r is specified at every grid point. Bicubic splines are used to find 
the interpolated value at any intermediate point. We can change the parameters 
of the bucket and see the effects on the torque developed. This anaysis will help 
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selecting the optimum shape of the bucket. This work was partly supported by a 
project sponsored by BHEL, BHOPAL. 
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Chapter 2 

Theoretical Formulation 


2.1 Introduction 

The water jet striking the Pelton Wheel buckets is modelled as a series of water 
paricles. In this chapter, the kinematics and dynamics of the fluid particle for tw r o 
dimensional (2D) and three dimensional (3D) situations are formulated. First we 
consider a 2D model for calculating the force and torque imparted on the rotating 
vane by a single particle. A general shape of the vane whose equation is given 
in polar coordinates is considered. Later, this 2D model is extended for a series 
of continuous particles and the concept of mass flow rate is introduced. Finally 
the general equations for 3D situations with any general shape of the bucket are 
obtained. 


2.2 Kinematics and dynamics of a single par tide for 
two-dimensional vane 

Let us consider a two-dimensional vane rotating with a constant angular velocity 
Lok about the point o as shown in figure 2.1 . A moving coordinate system having 
its origin at o’ with oo’=R attached to the vane. The equation of the vane profile 
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Figure 2.1: Schematic diagram for 2D case at the time cf impact 


in this moving coordinate system is described in polar form p = p(6). Initially, 
at t — t p o, let a particle (of water ) strikes the vane at the point with 6 = 6 0 . 
At this instant, the incoming absolute velocity of the particle is given by Vj. The 
trajectory of this particle before striking the vane is horizontal and at a distance 
H from the axis of rotation (o) of the vane. The inclination of the arm from the 
vertical is Vo ( figure 2.1 ) which is known. It is easy to calculate 6 0 from geometry 
as follows from figure (2.1), and A o'ap b0 


cos Vo _ cos (Vo + 00 ) 

TOW - -da-R ' 


(2.1) 


Equation (2.1) is solved numerically to obtain 6 0 from known values of H, R, 
Vo and the given functional form of p(9). 

For all subsequent instants, we assume that the particle follows the vane profile 
until it loses contact, when the force exerted by the particle falls off to zero or 
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becomes negative. To obtain the location of the particle on the vane at subsequent 
time instants, i.e. 6(t) with the known values of 6 0 , we need to solve particle’s 
dynamical equations as detailed below. 



Figure 2.2: Schematic diagram for 2D case after impact 


Referring to figure (2.2), let the point p b on the vane indicates the location of 
the particle at any instant t. The point c’ denotes the centre of curvature of the 
vane profile at the point p b . The absolute acceleration of the particle at this instant 
is given by 


Gabs = ao + u x (63 x p) + 2u; x (pp + p66) + (p — p8 2 )p+ ( pO + 2 p0)8 (2.2) 


where, the acceleration of the point o’ is 


7T 


ao = (jJ 2 Rcos(tt — d)p + uj 2 Rcos(— — 6)6 


(2.3) 
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Thus, 


&abs = (—oJ 2 Rcosd — io 2 p — 2ujpd + p — p9 2 )p+(u: 2 Rsind + 2up + p6 + 2dp)d 
= a p p + a e 0 (2.4) 

The components of a a b s along n and t directions ( the normal and tangent to 
the vane) are 


and 


a n = —a p sin <j> + ag cos 6 


a t = a p cos (j> + a$ sin <p 


where q , the angle between p and i is obtained from 


cos <f> 


+ 7 2 


(2.5) 

( 2 . 6 ) 


(2.7) 


p = pd 


( 2 . 8 ) 


and 

p = p'e + p'e 2 (2.9) 

where p' is the derivative of p w.r.t. 9, and 6 and 6 are the first and second 
derivative of 6 w.r.t. time (#). 

With the assumption of frictionless vane, it follows that the tangential accel- 
eration a t given by equation (2.6) is zero. Using equations (2.7}-(2.9), a t = 0 
reduces to 


6(-u! 2 R cos 6 + p' cos <b) - ( uj 2 R cos 8 + ui 2 p -I- 2 xp6 — p"9 2 + pd 2 ) cos <f>+ 
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(u 2 R sin 8 + 2up'8 -T 2 p8 2 ) sin o = 0 


( 2 . 10 ) 


Equation (2.10) is a second order differential equation in 8 only as q> can be ob- 
tained from the vane profile [see equation (2.7)] The ODE can be solved by using 
the fourth order Runge Kutte method (RK4) which will be discussed in detail in 
the next chapter. Since equation (2.10) is of second order, we need two initial 
conditions 0 O and 8 0 (i.e., the value of 0 and 6 at t = t pQ . 

The value of 8 o is known from equation (2.1) as stated earlier. To obtain the 
value of 8 0 , it is assumed that after striking the vane, the velocity of the particle 
along the direction normal to the vane (n ) is given by that of the point p b0 along 
the same direction, whereas the component of vj along the t direction remains 
unchanged. This assumes that the impact is perfectly plastic but frictionless. 

Thus at t = t p0 0 + (i.e., immediately after the impact ), the velocity of the 
particle (see figure 2.1 ) 


i% bs = uop b cos(0 o + 8 0 - A 0 )n + Vj sin(i 9 0 T 0 O + ipo)t 
= (coopb cos(<j>o + 8 0 - A 0 )(- sin p 0 ) -1- Vj sin (6 0 + <j> o) 

COS 4>q)p + ( UJOp b COs(l9o + <?0 - Ao) cos 00 + 

UjSin(0 o + + 0o) sin <£>o)^ (2.11) 


where A 0 is obtained as explained below. 

From A oo'pbo in figure (2.1) using sine laws we get 


R 


P o 


opb o 


( 2 . 12 ) 


sin(# 0 — A 0 ) sin A 0 sin 8 0 
From equation (2.12), given 8 0 , p 0 , R, we can solve for Aq and the distance 


OPb o- 

Further, the absolute velocity of the particle at any instant can be written as 
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V a bs — Vq F U! X p + pp + p09, 


(2.13) 


where the velocity of the point o’is 

Vo = ujRcos(— — 9)p + uRcos90. (2.14) 

Using equation (2.14)in equation (2.13) and substituting t = t p0 0 + , i.e., 9 — 

9q , 9 = 0o, we can write 

v° abs = (u)R sin 9 0 + p 0 )p + (uR cos 9 0 + wp 0 + p 0 9 0 )9. (2.15) 

Now, the (9-component of v% bs from equations ( 2.11 ) and (2.15) we get 

• {ojop b cos(0 o + <f) 0 - A 0 ) cos <p o + Vj sin(0 o + <?o 4- Tpo) sin 4>o - cupo — uR cos 9 0 ) 

9o = 

Po 

(2.16) 

which along with 9 0 gives the initial conditions for equation (2.10). 

Thus, after solving equation (2.10) numerically, we get 0, 9 and 9 at any instant 
t. Using these values, we get a abs from equation (2.4 ). In fact, only the normal 
component of a a f, s is sought ( &abs =0). The force exerted by the single particle is 
given by 


F n = —m p a™ bs (2.17) 

where, m p is the mass a single particle (= m p At, where rn p is the mass flow 
rate of the jet ). 

The initial force exerted by the particle during the first impact at p^ is obtained 
as 


F? = - %)" (2-18) 

After obtaining 0(t), we can track the particle in the fixed frame with the origin 

IUrraAL LIBRA#* 

». I.T., KANPUB 
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at o by writing its coordinates as 


y(t) = i?cos(^ 0 + ut) + pcos(6 + tp 0 + ut) (2.19) 

x(t) — R sin ("0o d - a>f) 4* psin($ + 4" cut) (2.20) 

2.2.1 Multi-particle system 

We consider the water jet as a series of particles. That is, particles keep striking 
and travelling along the vane one by one until they lose the contact from vane. 
At a given instant of time, a number of particles will be present on the vane and 
each particle has its own independent history. The particle forces are computed 



tpO 1 j *P e 


Figure 2.3: Multi-particle system 

independently. At each time interval, for the given position of the vane, a particle 
is allowed to impact the vane. The entire time-history of that particle is then 
computed, along with the rotation of the vane, till the particle loses contact with 
the vane. The next particle is now calculated, starting with vane’s pasition that is 
after the previous position. This particle’s time history is again calculated. 
Finally, after all the particles’ histories are computed, the results are collated to 
find the force, say, at a particular time t. 

Thus from the seperate particle computations we obtain F xpi , the force exerted 
by the particle at time U. After all the particle computations are done we know 
each F xpi for t p0 < U < t pe , where t p0 and t pe are times at which the particle p 
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impacts, and leaves, the vane respectively. To find the force F X i on the vane at 
time ti, we have to find the sum 

F xi = ^2 Fpi f or °-tt P suc h that t p0 <U< t pe 

p 

A similar approach can be taken for the torque. 

2.2.2 Determination of initial position of a paritcle 

As shown in figure 2.4, the first particle strikes the vane when the arm makes an 
angle ip 0 from the vertical and keep on striking one by one until the vane rotates 
away from the line of the jet. The vane rotates with angular velocity oj and angle 
ip made by the arm with the vertical at any instant of time is ip 0 + cut. The last 
particle strikes the vane when the arm has rotated to ip max given by 



Figure 2.4: Initial and final position 
ipmax = arccos( cos( 0 o ) ) 

For a given position of the vane, ip, we want to find the 6 where the particle 
will strike. From this, initial position of the particle path will be computed. From 
figure ( 2 . 2 ), in A o'api 

cos ip _ cos (ip + 9) 

= 2fF* 
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or, 

p{ 9 ) cos(ip + 9 ) + Rcosip — H = 0 (2.21) 

With given p( 9 ), the value of 9 satisfying equation ( 2 . 21 ) with ip = + ut 

gives the location at which the particle strikes. 

2.2.3 Particle torque at the time of impact and after the impact 

Momentum is imparted by the particle when it strikes the vane and also during the 
subsequent motion until the particle loses contact with the vane. Initial impact for 
a single particle having mass m p (= m p At), where m p is mass flow rate of jet , is 
given by 

f initial = ^pi^abs ) 

= -m p (ojop b cos(6> 0 + (p Q - X) - Vj sin(# 0 + <p o - f ) 

So torque imparted by the particle about the point o at the time of impact 

M initial = ^initialOPb Sm 70 

= -m p (u>opl sin7o(cos ( 9 0 + <j ) 0 - A) - Vj sin(6 > 0 + <j>o ~ |) 
with 7o = ^>0 + - A 0 - f 

The torque imparted by the particle about the point o at any instant during the 
motion of the particle on the vane is given by 
= —m p a n opb sin 7 

To calculate the total torque imparted by the particle about the point o during 
its motion by the vane, we track the particle at every instant of time and sum up 
the torque to yield 

£ M At = -m p 12 a n op b sin 7 

We can obtain 7 and opb from the A 00' p b in figure 2 . 2 . Using sine laws we 
get 

' R = 0.22) 
sin (0 - A) sin A sin# 

and 7 = ^ + # — A — | 
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So the net torque imparted by the particle is summation of the initial torque 
and on subsequent motion. 

Thus, the net torque= £ M Af + M initial 


2.3 kinematices and dynamics of a single particle 
for 3-dimensional bucket 


Referring to figure (2.5) a, let us consider a bucket of general shape rotating with 
a constant angular velocity uj about the X-axis. There are two frames of refer- 
ence considered, one fixed frame XYZ centered at o while another one moving 
frame xyz centered at o’ .The frame xyz is attached to an extended point (o’) of 
the arm whose distance from o is R. In figure 2.5 (b), the sign convention for 
spherical coordinates is shown. The surface equation of the bucket is assumed to 
be ip(r, 6, <p) = 0 and when its gradient is given by 


7* dip „ 1 dip p: 

Vv= d^ r + rm e + 


1 dv : 
r sin 8 do ' 


Initially at t = 0, the arm is assumed to make an angle X'o from the Z-axis, 
a particle travelling along the Y-axis strikes the bucket at point pw with velocity 
fij.The postion vector of p w with respect to the moving frame is rf. After striking 
the bucket at point p w , the particle is assumed to travel along the bucket. At 
that instant, absolute acceleration of the particle is along —Vip (i.e. the bucket is 
assumed to be frictionless so the particle is subjected only to a normal force and 
expreriences a normal acceleration) and given by 


a abs = o 0 -l- u x (£ x r) + 2lj x (rr + rdd + r sin 8<p6) 4- a T f + -1- a ^ (2.23) 

where Go is the acceleration of point o’ 

Gq = —uj 2 Rk 
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and a T , a$, a# are the acceleration components of the point Pb in f, 8, o direc- 
tions, respectively. These components are given as 
a T = r — rd 2 — rcj) 2 sin 2 8 
a$ = rd + 2 rd — rq 2 sin 8 cos 8 
a# = rep sin 8 + 2 rep sin 9 + 2r8<p cos 8 

The relation between the cartesian coordinates and spherical coordinates for 
the moving frame are given as follows 
f = sin 8 cos (pi + sin 8 sin (pj + cos k 
8 = cos 8 cos (pi + cos 8 sin (pj — sin 8k 
4> = — sin (pi + cos (pj 
The above relation can also be written as 


< 


1 

sin 0 cos (p 

cos 8 cos <p 

— sin (p 


r \ 

f 

= 

sin 8 sin (p 

cos 8 sin 4> 

cos<p 

< 

8 > 

1 

cos 8 

-sin# 

0 


4> 

\ / 


So we finally write equation (2.23) as 


Q-abs — A^r -f- Ag9 "T AqQ 


where. 


(2.24) 


(2.25) 


A T 


—u) 2 Rcos6 — cu 2 r(cos 2 #cos 2 <p + sin 2 <j>) — 2u;r(<£cos#sin# 
cos (p + 9 sin (p) + (f — rd 2 - rep 2 sin 2 9 ), (2.26) 


Ag 


oj 2 Rsmd + u 2 r sin 8 cos 8 cos 2 <p + Cor sin0 + 2u>(r<^sin 2 # 

+r sin 4>) + (r# + 2f# - r(p 2 sin 0 cos 9 ), (2.27) 


and 


17 



Acj, = — w 2 rsin 9 cos (f>sln <fi + uir cos 9 cos (j> + 2w(r cos 0 cos 4> — 

—rd sin 9 cos p) + (r <j> sin 9 + 2r ^ sin 6 4- 2r9o cos 9) . (2.28) 

The following two equations are obtained by assuming the bucket to be fric- 
tionless. 


and 


along with 


Ar 

dip 

dr 


Ae 

1 dip 
r dd 


(2.29) 




A r 

dip 1 dip 

dr r sin B d<p 


(230) 


4>{r,6,<f>) = 0 (2.31) 

Substituting equation (2.31) in equations (2.29) and (2.30) we get 


flu 9 + Qi2$ + ^13^ 2 + CL\a4^ "T Ul5<^ "T CL\^9^) + C = 0 


(2.32) 


and 


021^ T- 0-22.0 "T &2z9~ + + d25& + 0-26^^ + d — 0 (2.33) 


where, 
flu 


(_l(~e D 2 

V r 


di* d\ 

s* l do do 

a l2 — ~r~^~ 

_ f dv \2 <9^ V 

Q — g [MpL&x&t 
13 r (&)* 


r 5v' 

1 dr ■ 


I .fit 


d‘ v r d%_ >3 
Tv dO > 


(£) 3 


d^t/> djb 

i mZJL _i_ 

r 2£ h 
dr 

t»2,. 


as 


an 


d^ v dv dt d^v dv ( dtp \2 ^ 'PLQjsL 

l - 1 g? % f ] - ll&g- - sin 2 e d 4 A r sin 9 


<W 


(f) 3 


dv 

757 
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+r sin 9 cos 8^ 


ai5 = — 2a; cos 0 sin# cos 3 §£ - 2o>rsin 2 #cos3^|: + 2o;sin3§§ 


qLh ( Oil gPtI> 5t!» 


Ol6 


2 3r36 36 33 , 2 Qrdc ^ 36 ) 

t dtjj\ 2 ‘ r ' ^ 3i p 


+2g 


- ^_Z_ 

2 3033 30 

^ — 

a7 


dr 

3 2 3 / 33 \2 3v 
2 d73~ v de ) do 

~ r (§7p 


c = — pff (u) 2 Rcos9 + a’ 2 r (cos 2 # cos 2 3 + sin 2 3)) - |^(a; 2 i?sin0 
+u 2 r sin 9 cos 9 cos 2 6 + Ccr sin 3) 

021 = 


33 3v 
de dc 


" r sin 3 ST 

OT 


022 


1 


(§£P 


- rsine - ~^ -- rSing | : 


_ 2 

0,23 ~ rsin 0 


3 2 3 30 3v 


>r30 36 3c 


1 


TUP 


3 2 tw 3Q , %2 3^ 

de > 33 


(¥p 


33 _ 


3 2 t/> 33 


0-24 


jy^flP 

lfP~ 


rsin 3 
3 2 3 33 


1 dll' 
sin & d(p 


3-v / 3y >3 x- ret. _ , 

1 3t- v 3o l 1 33^ <93 1 _• n dv 

r sin 6 rsin 3 


rsin 3 (^)2 rsin 3 (§-) 3 rsin# 1 ,JAAJ ' v d$ 

a 25 = + 2a-’ cos 6 cos 3|f + 2a;r sin 6 cos 3§f 

— 2 l^glffP , 

26 r sin ® (^P + ■ — 

^H(HP 


ae 

3 2 v 3Q‘ 3Q> 
2 3roo 36 3 q 
r sin 6 (|£)2 


rsin 3 


3 2 3 3t/' 
3633 c?3 


33 

3r 


rsin 3 


— - 2r cos 9§ + 2 sin 9% 
d = uj 2 R cos 9 - a; 2 r( cos 2 # cos 2 3 + sin 2 3)) - §J:(— w 2 r 

sin 6 cos 3 sin 3 + wr cos 6 cos 3) 

Equation (2.32) and (2.33) are two coupled second order differentia] equations 
in 9 and <f>. So we need initial conditions 9 0 , 3o> 0o and 3 0 (i.e. the values of 9,9, 
3 and 3 at t — f p0 ) as discussed in section 2.2. The value of 9 0 and 3o can be 
calculated from the following equations. 


X c = 0 

l' c = R sin xo + 2 sin Xo + V cos xo 


Z c = H = R cos xo + ^cosxo - V sin Xo 
where xo is the angle made by the arm from the Z-axis when the particle 
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strikes. The methodology to solve the above equations and to get # 0 and <fr 0 will 
be discussed in detail later on. 

To find the value of # and 4> at t = 0, we proceed as follows: 

As discussed in the previous section, after the impact the component of jet 
velocity Vj along Vtp is same as the velocity of the point p b in that direction but 
the velocity of the jet is unchanged in the direction normal to Vip, i.e. tangential 
to the surface. 

So after the initial impact at t = 0 + , 

along Vip, the component of Vj =component of v pb 

or. 



UJ A 

W 


Rk + oj x rf + rf -f r99 + r sin 9(j)<f) 


(2.34) 


where 

Vj = Vj cos xj + Vj sin \k 

The relations between the cartesian coordinates of the fixed frame and the 
moving frame are as follows; 

} = % 

J = cos xj + sin xk 
K — — sin xj + cos xk 

So from equation (2.34) we equate components on the left and right to find 


— ojR cos 9 sin 6 
0 O = 


and 


jv^Frfe + (cos x cos# sin d - sinxsin#) — wrsin# 

(2.35) 


4>o 


-ioR cos 4>-ur cos 9cos4> + Vj cos % cos <p + 


drl’ 


|V^| 2 rsin 9 d<p 


rsin# 


(2.36) 


20 



where, 

Vj .V% p = Vj( (cos x sin 0 sin $+sin x cos 6) ~ + (cos x cos 9 sin <j >— sin x sin 9)~^g 
+ cos x cos 4> — \-s ) 

Using these initial conditions we solve equations (2.32)and (2.33) numerically 
we get 0, 9 ,0 ,4> ,<j) and <j> from the initial (impact ) time t pi to the time the particle 
leaves the bucket. The particle is considered to have left the vane whenever the 
computed a n becomes positive. Using these values, we get a abs from the equation 
(2.25). The force exerted by the single particle is given by 

F = -m p a abs (2.37) 

where, m p (= m p At ) is the mass of the particle. 

The initial impulse exerted by the particle during the first impact at p b0 is 
obtained as 


I n — ^p( 



and the average force is taken to be 


m\ 


where 

v pb = u x op b , 
opb = Rk + rf 


Fn — 17lp{ 


Vj.Vlp 

w 


|V^| 


(2.38) 


(2.39) 


2.3.1 Initial positions for 3D case and modelling of water jet 
shape 

In section 2.2.1 the multiparticle system has been introduced and for 2D case. S- 
ince particles strike one by one and keep on sliding along the vane till they lose 
contact from it. The computation for each particle is done seperately, from differ- 
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ent initial conditions. For a particle on centre line of the jet, the initial position 
where the particle strikes is: 


X c = 0 

Y c = R sin x + z sin \ + y cos x 


Z c = H = R cos x + 2 cos x - 2 / sin x 

where, (X c , Y c , Z c ) is the point of impact with respect to fixed frame and 
X=(xo + ut) is the angle shown in figure ( 2.5 )which varies from Xo to Xmax, 
and (x c , y c , z c ) are the coordinates of the point of impact on the bucket surface in 
the moving coordinates system (figure 2.5 ). 

It is however unlikely that particle will come only along the centre line, as the 
jet some finite radius l 0 . So we assume that the particle may come along a point 
randomly located on the circular cross-section of the jet. We do this by assigning 
a random position to the point by the prescription: 

10 = radius of water jet 
b=jradius from o" 


Figure 2.6: water jet 



l = l 0 \/ranl 
ft = 2~ran2 

where rani and ran 2 are two different uniformaly distributed random num- 
bers having values between [0,1] , the random point at any l and ft are shown in 
figure ( 2.6 ). The square root in the prescription of l is to ensure that the points 
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are randomly distributed on the cross section.The initial position of the particles 
will be perturbed from the centre line by : 

AX = l cos /3 
A Z = l sin 

where AX and A Z are the variation in X-direction and Z-direction. respec- 
tively, from the centre line. So now the equations for the initial positions are 

X p = AX 


Y p = Rsin x + 2 sin x + y cos x 

Z p = H + AZ = R cos x + z cos x ~ V sin x 

where (x, y, z) is the initial point in the moving frame which has to be deter- 
mined for a given AX, A Z, %,H and R. This requires the use of of a root finding 
algorithm which will be explained later. 

2.3.2 Particle torque computation in 3D case 

In the same way as discussed in section 2.3.2 for the 2D case, we can calculate 
initial torque and on subsequent motion as follows: 

Torque imparted by the particle about the X-axis at the time of impact 

M initial = OPb X 

Torque imparted by the particle about the X-aixs during the motion on vane 
M& t = m p (op b x tip) 

= mp^RsmdA^f — (Rcos9 r)A,p9 + (Ag{R cos 6 + r) + 

yf r i?sin 9)(j)) 

£ M& t = £ rrip{op b x a p ) 

^p£(0P6 x d p ) 
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So net torque= £ M At + M imtwl 

The total torque contribution of all particles can also be obtained in the same 
way as discussed in section 2.2.1. 
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Chapter 3 

Algorithm and test cases 


3.1 Introduction 

This chapter presents the algorithm for particle tracking, force & moment calcu- 
lations and for the integration of many particles information to obtain the overall 
torque. Towards the systematic development of the software, first of all the algo- 
rithm is developed for the 2D case and subsequently for the 3D case. 

For each particle the path of the particle is computed by use of Newtonian 
dynamics. Subroutines were developed to compute the force exerted by the parti- 
cle and its moment of momentum at time of impact and subsequently during the 
motion on vane. The torque exerted by all particles is summed up and related 
to the mass flow rate. The computer program is used to simulate problems with 
known analytical solutions and results were found to be perfectly matching for 2D 
circular and non-circular vanes. 

The computer program was generalised to the 3D problem. The computer 
program was used to simulate the analytical results obtained previously and a 
perfect match was again found. The subroutines to compute torque etc., were 
extended to the 3D case so that the program can calculate, given the bucket shape, 
mass flow rate, the torqe developed by the Pelton wheel. 
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3.2 Runga-kutta (R-K) methods 


All the differential equations derived in chapter 2 are solved by using the fourth 
order Runga-Kutta 4 (RK-4) method. So RK-4 methods need a brief mention. 

Runga-Kutta (R-K) methods are a one-step procedures which require only e- 
valuation of first order derivatives, but which produce results equivalent in accu- 
racy to the higher-order Taylor formulas. These methods are used to numerically 
solve ordinary diffemtial equations (ODES) of first order. Any ODE of higher 
order n is first reduced to a set of n first order ODE’s and then solved. 

A first-order initial value problem (IVP) is of the form 


dy 

dx 


= f(x, y)\ y(x 0 ) = y 0 ] x 0 < x < x r 


(3.1) 


Since, in general, the analytical solution cannot be obtained, the true solution 
y(x) is approximated at N + 1 evenly spaced values of x, (ar 0 , x\ , • • • , x^), so that 
'h' the step size is given by 


h = 


x 0 


N 


(3.2) 


Xj-j-i — X{ "4" h (3*3) 

where i = 0, 1 , - - * , AT 

Beginning with the initial value of yo at the values of the solution 2/i , 2 / 2 , - - - at 
the points X\, X 2 , ■ • ■ are computed sequentialy according to the following rule: 

Vi+i = : yi + + 2^2 + 2A' 3 + A' 4 ) (3.4) 

b 

where 

Ki = f(xi,yi ) 


K 2 = ffa + ^h, yi + | hK x ) 
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•^3 — f{ x i + -h, yi + -fiJT 2 ) 

K 4 = ffa + h,yi + hK' 6 ) 

Consider the following system of n simultaneous first order ODE’s: 


dy i 
dx 
dy 2 
dx 


fl (^) ?/l j Vli ’ i Vn) 

hi*, Vl,V 2 ,’--,Vn) 


dy n 
dx 


fn (x, yi, y 2) - ••,?/«) 


with initial conditions 


( 3 . 5 ) 


yi(*o) = yi,o 

^2(^0) = ^2,0 

Vn(Xo) = 1/n.O 


The solution of such a system is along the same lines as the solution of a single 
first order equation repeated n times. The R-K algorithm is applied to each of n 
equations in parallel at each step. 

Consider a single higher order ODE 


d m y dy cPy d™ ^ 

dx m ~ F ( X ' y, dx’dx 2 ’'''' dx™- 1 * 


( 3 . 7 ) 
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with approriate initial conditions for 

dy[x o) d m ~ l y{x 0 ) 

''W- 

The above ODE can be rewritten as a system of n simultaneous first order 
differential equations. Then the RK4 algorithm can be applied as above. 

3.3 Two-dimensional simulations 

Figure 3.1 shows the simplified schemetic diagram of 2D case, having circular 
radius, at time of impact and after the impact. 



Figure 3.1: schematic diagram for 2D simple case 


As discussed in section 2.2, we have a single second order ODE ( after the 
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assumption that the vane is frictionless ) which is given as 
6(-oj 2 R cos 9 + p' cos <j>) - (u 2 R cos 9 + oj 2 p + 2up9 - p"9 2 + p9 2 ) cos (j>+ 

(u 2 R sin 9 + 2top'9 + 2p9 2 ) sin ^ = 0 (3.8) 

To track the particle at every instant of time, we need two inital conditions 9 0 - 
and 9 0 . The initial 9 0 is calculated as discussed in section 2.2. With the assump- 
tion, as discussed in section 2.2, that the velocity of the particle along the normal 
direction to the vane is given by that of the point p b0 , whereas the component 
along the tangential direction remains unchanged from the jet velocity, we get 

■ (uiopb cos(9 0 + (f>o — A 0 ) cos fa + Vj sin(# 0 + <f> o) sin (f> 0 - up 0 - uR cos 9 0 ) 

9q — 

P o 

(3.9) 

So the computation is done in following way: 

1. compute the inital conditions ( 9 0 ,9 0 )for the particle. 

2. compute the initial impact force imparted by the particle. 

3. compute the value of 9 and 9 at next time step using previous step’s 9 and 9 
and the RK4 solver. 

4. compute the force and torque by the particle using the computed 9 and 9. 

5. step 3 and step 4 repeated until the force exerted by the particle falls off to 
zero. 

i 

3.4 Test case for RK4 solver 

We assume the simplest case when there is no vane present to obstruct i.e. the 
particle is like a free particle as no force will be exerted on the particle, its absolute 
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acceleration will be zero and particle will move in a straight line. 

In our numerical simulation we suppose a fictitious vane (see figure 3.1) hav- 
ing circular vane of radius p - c. So we solve equation (2.10) under following 
two initial conditions: 6 = 8 0 and, from equation (2.16), 9 0 = 
since ip, 0 = 0 and <pu-\. 

For different values of c, 0 O . Vj, u and R the equation (2.10) is solved having 
the time step At = .01. After solving the equation (2.10) for 8 we track the parti- 
cle at every' instant of time using equations (2.19) and (2.20) and the path travelled 
by the the particle was found to be straight line. This validates the kinetics of the 
panicle. 


3.5 Test case2 

This problem is selected to test the numerical solution against an analytical solu- 
tion in a complex dynamic problem. 



Figure 3.2: Schematic diagram for test case 2 


A frictionless circular ring having diameter 2c and centre at o is hinged at 
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point o. Initially a particle is at point p b0 at the diametrical distance from o as 
shown in figure ( 3.2 ) and the ring is at rest. At t = 0, the ring starts rotating 
counter clockwise at constant angular velocity co 0 . At any instant of time the 
particle is at point p b as shown in figure ( 3.2 ) whose position (9), angular velocity 
(6) and radial force (a p ) can be obtained analytically and are given as follows: 


9 = -4arctan(exp(w 0 t) - 


7T 

4 


) 


9 = -w 0 \/2(l + cos#) 


dp — ^qC 


2y2(l 4- cos 9) — 3(1 + cos 9) 


The situation shown in figure 3.1 can be reduced to the situation as described 
in figure 3.2 by applying the following conditions: Vj = 0, u = -ui 0 , p = c and 
R — p = 0. So equation (2.10) is solved using RK4 solver with following initial 
conditions: 

0 O = 0 and 9 0 = 0 

An angular velocity ui=2, c=5 and time step At = .01 are used for computa- 
tions.The results are shown in the following graphs. 

The following graphs shows the perfect match between computed solution and 
analytical solution for 9 , 9 and a p . So this test validates the numerical solver for 
the equation ( 2.10 ) completely. 


! 

i 
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3.6 Test case 3 

This problem is selected to check the computed force and torque applied by a jet 
against known analytical solutions for a 2D vane. An anlytical solutions are not 
known for even the simplest cases of a rotating vane, we have considered two 
special cases for which such solutions are known: (a) a stationary vane and (b) 
a translating vane. These simple cases however allow us to check the numerical 
solution which is computed by essentially the same method (solving the ODE’s 
etc )as are more realistic cases. The utility of this test is that in this case we check 
not only the particle dynamics but also the integration routines which collate all 
the particles information to give a composite picture of the force and torque on 
the entire vane. 

Referring to figure (3.6), the vane has the shape of p(6) = 1+cos 9. The water 
particle strikes the vane at point ’a’ and after travelling along the vane it exits at 
’b\ 

case a: stationary vane 


33 





o 


Figure 3.6: Schematic diagram for test case 3 

For this case, the vane remains in the position shown in figure (3.6). The 
analytical solution for the force exerted along the X-direction ( F x ) and the torque 
about the point ’o’ are easily calculated as 

F x = rhvj( 1 + ^=), 

and 

T = mvj (^~-) + (R + 2) . 

The program is simulated with following conditions: u> = 0, m = 1 and initial 
conditons 0 0 = 0, p = 2 and 6 = ^ sin ^ 0 ~-— • Computations are done for 
different value of R and Vj and are given in the table (3.1). 

Table (3.1) shows the perfect match between the analytical and computed solu- 
tions. The time step At is very important, for if particle comes with high velocity, 
it is going to take very little time to travel along the vane. So time step has to be 
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velocity(m/s) 

R(m) 

force(N) 

analytical 

force(N) 

computed 

torque(Nm) 

analytical 

torque(Nm) 

computed 

At 

time step 

1 

1 

1.707200 

1.70836 

3 

2.99881864 

6.01 

5 

1 

. — - - 

8.536000 

8.53928 

15 

14.9941990 

0.01 

50 


85.36000 

85.42057 

320.721344 

320.791016 

0.001 

100 

" 3 

170.7200 

170.8534 

641.628784 

641.442688 

0.0001 


Table 3.1: For stationary vane 


smaller than the time step taken for relatively slower velocities, 

case b: moving vane 

When vane in Fig. (3.6) has a translational velocity v 0 . Then the force F x is 
given by 

F x = t h p (vj - Vo ) (1 + 

The computer program is simulated using following conditions: 

R — > oo, u) — > 0 and intial conditions 0 O = 0, po = 2 and 6 — V] sm v ° 
and computation is done for different value of Vj, v 0 and m p and are shown in the 
table (3.2).Thc numerical solutions match very well with the analytical solution 
in this case also. 


Vj (m/s) 

v 0 (m/s) 

rh p (kg/sec) 

fofce(N) 

analytical 

force(N) 

computed 

At 

2 

1 

1 

1.707200 

1.7078361 

0.01 

10 

5 

5 

42.68000 

42.71034 

0.01 

50 

20 

5 

256.0800 

256.5815 

0.001 

100 

30 

c 

-J 

597.5200 

599.8303 

0.0001 


Table 3.2: For moving vane 


The equations ( 2.32 )and ( 2.33 ) derived for 3D situations were also simulated 
for the above test cases and a perfect match was found with the 2D numerical 
solution as well as with the analytical solution. 
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3.7 Interpolation scheme for the suface of the buck- 
et 

The shape of a Pel ton wheel bucket is quite complex and cannot be represented 
by simple a analytical function. To overcome this difficulty, the bucket surface 
is represented in the form r = f(9, 4>) wher data are supplied in the form of a 
rectangular grid ( M x N) in 8 and 4> and at each grid point the radial distance (r) 
is given. To calculate the value of r at any 9 and (j) we need some interpolation 
scheme. 

We use the cubic spline interpolation scheme since it provides a superior ap- 
proximation for functions that may have local abrupt changes. The objective in 
cubic splines is to derive a third order polynomial for each interval between grid 
points. In one-dimension the function would be 


fj(x) = QjX 3 + bjX 2 4- CjX + dj (3.10) 

where j(= 1, 2, 3 •••,«) 

To develope cubic splines we have to solve 4n equations for the values of Oj, 
bj, Cj, dj for n intervals. These equations are obtained by imposing the condi- 
tions that the function yeild the tabulated values at the grid points and also that 
the piece-wise functions have continuous first and second derivatives at the grid 
points. 

The second derivatives of the above function can be represented using a first- 
order Lagrange interpolating polynomial as 


/: (x) = f (Xj.r) 


X Xj ^ pH 
Xj—i Xj 


(Xj) 


X Xj — j 
Xj Xj — i 


(3.11) 


where /" (x) is the value of the second derivative at any point x within the j th 
interval. Integrating the equation (3.11) and using function equality conditions 
i.e. that f{xj- 1 ) must equal the tabulated fj and f(xj-i) must equal /j-i,and 
invoking the conditions that the first derivatives at grid points must be equal i.e. 
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= f'n iM.we get 


/ (*, -t 1 ) -f (£V) /fa) - /fa-i) 




x j+l 


Xi 


Xi 


Xj-i 


(3.12) 


These n-2 linear equations have n unkowns, all the f“(xj)'s. So setting f"(x i) 
and f"(: x n ) to zero, which gives the so-called natural spline which has zero sec- 
ond derivatives its boundaries, the above equations can easily be solved by using 
tridiagonal algorithm (TDMA). 

In two dimensional interpolation, we seek an estimate of f(9 , <ft) from 2- 
dimensional rectangular grid (M x N ) of tabulated values of each of each of 
the independent variables 9, <f>. To interpolate one functional value , i.e., of r, we 
perform M one dimensional splines evaluation across the rows of the table, fol- 
lowed by one additional one-dimensional splines down the newly created column. 
This is called bi-cubic spline interpolation. 

Our primitive equations (2.32) and (2.33) for 3D computation contains the 
terns fjf, g£, 0 and g|. As discussed earlier we have reduced 

V(r, e, *) lo r - /(«, *). So g = 1, 0 = 0, = 0, ff = % = -g, 

00* - Qpi - ana d6d(j> dBdf 

To obtain the first, second and cross derivative with respect to 6 and <f>, at any 
point we apply finite difference scheme as we know the value of f(9,<l>) at any 
point on given domain by using bi-cubic spline interpolation. So we can write 

86 ~ 2A 6 

anri hMi h zMAi )- etc. 

auu gp — A0 2 

where A0 is a very small value. We use bi-cubic splines to evalutate all the 
functional values f(9,<j>), f(9 + A 9,4) etc and finally get the first and second 
derivative with respect to 9. Similarly we can calculate first and second deriva- 
tive with respect to <j> and also the cross-derivative. The computed value of first 
derivative and second derivative of f{9, <f>) with respect to 9 and <f> is found to 
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be very close to the analytical solution for any smooth function but cross deriva- 
tive (&&) cloes not match wcl1 - This is a fundamental problem with the bicubic 
spline because this icpiesentation does not account for cross terms. Since in e- 
quation (2.32) and (2.33) cross derivative comes only once as summation term, it 
doesnot affect the final solution very much. However, it is reccomended that the 
program should be improved by using rational B-splines to represent the bucket 
surface instead of bi-cubic splines. 

3.8 Calculation of initial postion for the series of par- 
ticles 

Before the calculation for a particle can start we need to know its initial condition- 
s, i.e., its initial position and velocity, both radial and angular need to be known. 
For a vane of complex geometry the first difficulty that arises is to find where the 
particle strikes the vane. As the vane is constantly moving, it has a different posi- 
tion for each particle. To calculate the initial position of particles in the spherical 
coordinates of the moving frame, the following methodology is used: 

Referring to figure 2.5 (a) and 2.5 (b) the postion of the particle in the fixed 
frame of reference is given below as discussed in section 2.3.1. 

X p = AX (3.13) 

Y p = J? sin x + z sin x + y cos x (3.14) 

Z p = H + AZ = .R cos % -f 2 cos x ~ y sin X (3-15) 

We know the value of A A', H + AZ, x and R for every particle as discussed 
in section 2.3.1. Equations (3.13), (3.14) and (3.15) can be cast as follows: 


38 



(3.16) 


* x = X p 

y = Y p cos % - Z p sin \ ( 3 - 17 ) 

z = Ip sin x + Z p cos x ~ R (3.18) 

The equation of surface ip(r, 9,<f>) = 0 can be cast as ^(x, y, z ) = 0. 

Using equations (3.16), (3.17) and (3.18) we can cast the equation ip{x, y, z) 
0 as 

ip(X p , Y p , Z p ) = 0 (3- 19 ) 

In equation ( 3.19) X p and Z v are known, hence Y r is the root of the equation 
(3.19). It can be solved by using NEWTON-RAPHSON scheme. After getting Y p 
we can get 9 and <f> by back substitution. 

For 3D computation, we have represented ip(r,6,<j>) as r - anc * ^ ata 

are supplied in the form of rectangular grid (M x N) in 6 and <f> and at each grid 
point the radial distance (r) is given. To find the value of 8 and (j) at the time of 
impact for each particle we follow the follwing procedure. 

First of all, we convert the rectangular grid in 6 and (j) into the rectangular grid 
in X p and Z v having the value of Y p at each grid point using equations (3.13)- 
(3.15). For a given value of X p and Z p we interpolate the value of Y v using bi- 
cubic spline. Now we can calculate the value of x, y, z in moving frame using 
equations (3.13), (3.14) and (3.15)and from known value of x, y, z we calculate 
the value of r, 6 and d > . 

3.9 Discussion on a complete test case 

The following steps are taken for entire comptutation. 

step 1 - grid generation in 9 and (j) from given x, y,z points. 
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step 2 - computation of a number (n) of particles striking the bucket after 
assuming the time step ( At ) and starting angle of rotation (ipo ) for given angular 
velocity of bucket (u ). 

step 3 - initial condition (# 0 ,</>oA and 0 O ) of the particle is computed, 
step 4 - initial momentum is computed for the particle, 
step 5 - use the RK4 solver to compute the value of 9,<j>,0 and 4> for next time 
step. 

step 6 - force by the particle is computed and checked whether it is less than 
zero or not. If the force is zero that means it has left the bucket. 

step 7 - repeat step 5 and step 6 for the particles present in the bucket at that 
time step. 

step 8 - repeat step 3 to 7 till all particles strike the bucket, 
step 9 - collate all particle information to obtain the overall force, torque, pow- 
er, efficiency, etc. 

We have taken a test case for a hemispherical bucket. The parameters are as 
follows: 

velocity of jet striking the bucket Vj =10m/s. 

angular velocity of bucket w =3 rad/sec 

The height at which jet strikes H= 1.7m 

Time step At =.01 sec. 

mass flow rate m p =2 kg/s 

The bucket statrs from ipo=— f to | 

Number of particle striking the bucket n=34 
power develped by wheel =35.58 w 
Efficiency of the bucket obtained is =37.018% 

We tracked the every particle at every instant of time and torque is computed 
using above algorithm. The figure (3.7) shows time history of all particles hitting 
the bucket. Paticles strike the bucket at every time instant and after travelling along 
the bucket they lose the contact. So figure (3.7) shows that maximum number of 
particles present in bucket is no more than 25 at a time while number of particle 
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striking the bucket are 34. 


25 


O 

o oo 


s 20 

OQ 


w 


o o < 

o o o 

o oooo 


15 b 


CL. 

U3 

b 

F 

P£ 

2 

§ 

g 


10 


_1 L_ 


10 15 20 25 30 35 40 45 

TIME STEP 


Figure 3.7: Particles’ time history 


Figure (3.8) shows the torque history of a particular particle over entire time 
span, starting from the initial impact till it loses the contact form the bucket. Al- 
though the efficiency of the hemispherical bucket is very low but this test case 
shows the general behaviour ( as shown in figure 3.8 ) of the Pel ton wheel that the 
maximum energy is transferee! during the impact and slowly it falls off to zero. 
As discussed in chapter 1, the maximum energy can be transferred when the par- 
ticle exactly bounces back in the opposite direction. But to avoid the interference 
from incoming flow, particle is deflected through some angle. Here the surface 
geometry comes in the picture. 


Figure (3.9) shows the torque history of particle during its motion on vane. It 
clearly shows that a substantial amount of energy is transferred during the parti- 
cle’s motion on the vane. So depending upon the bucket shape, this energy can be 
enhanced. So a lot of improvement can be predicted on the basis of this analysis. 
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TORQUE IMPARTED ON THE BUCKET AFTER IMPACT 



Fig' 





3.10 Conclusion and future work 


In this thesis, the torque is calculated for a single bucket having simple geomet- 
rical configuration. The scheme can be used to test the complex geometry. The 
scheme can also be easily extended for multi -bucket system. 

In the future work, attention can be paid upon complexity of the geometry of 
the bucket like, the shaip edge in the centre of bucket and cut on the top of the 
bucket. The present work can be extended for non-friction surface. Attention can 
also be paid upon the water jet which is assumed to be running in straight line 
while in reality it has some curvature. 
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